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ABSTRACT 

We report results from calculations investigating stationary magnetic field configurations in accretion discs around magnetised neutron 
stars. Our strategy is to start with a very simple model and then progressively improve it providing complementary insight into results 
obtained with large numerical simulations. In our first model, presented here, we work in the kinetic approximation and consider the 
stellar magnetic field as being a dipole aligned with the stellar rotation axis and perpendicular to the disc plane, while the flow in 
the disc is taken to be steady and axisymmetric. The behaviour in the radial direction is then independent of that in the azimuthal 
direction. We investigate the distortion of the field caused by interaction with the disc matter, working in full 2D and solving the 
induction equation numerically. The influence of turbulent diffusivity and fluid velocity on the poloidal field configuration is analysed, 
including discussion of outflows from the top and bottom of the disc. We find that the distortions increase with increasing magnetic 
Reynolds number (calculated using the radial velocity). However, a single global parameter does not give an adequate description 
in different parts of the disc and we use instead a 'magnetic distortion function' Din(r, 6) (a magnetic Reynolds number defined 
locally). Where £)„, <c 1 (near to the inner edge of the disc) there is little distortion, but where » 1 (most of the rest of the 
disc), there is considerable distortion and the field becomes weaker than the dipole would have been. Between these two regions, 
there is a transition zone where the field is amplified and can have a local minimum and maximum. The location of this zone depends 
sensitively on the diffusivity. The results depend very little on the boundary conditions at the top of the disc. 

Key words. Accretion, accretion disks - Magnetic fields - Magnetohydrodynamics (MHD) - Turbulence - Methods: numerical - 
X-rays: binaries 



1. Introduction 



I Magnetic fields play a fundamental role in the physics of accre- 
tion discs. First of all they are thought to be the origin of the 
turbulence which makes the accretion itself possible: this turbu- 
lence is usually thought to be caused by the magneto-rotational 
instabihty MRI (Velikov il959l Chandrasekhai- 1960 Balbus & 
Hawley 1991) although recently other important instabilities 
have been suggested that can operate even when MRI cannot 
(sheer-driven instability, Bonanno & Urpin l2006l 120071 current- 
driven Tayler instability, Tayler 119731 Riidiger et al. 120071) . 
Secondly they can also determine the geometric and kinetic 
structure of the disc. In addition, magnetic fields are invoked to 
I explain several characteristic features of accreting systems, such 
I as particle coUimation (jets), radiation collimation (pulsar light- 
house eff'ect) and spectral line production (cyclotron and syn- 
chrotron emission). There are hints that magnetic fields are ac- 
tive even in the dead zones of protoplanetary discs where they 
transfer angular momentum outwards (Turner & Sano 2008). 
Overall, it seems that regardless of the particular kind of ac- 
creting system, from ones around supermassive black holes in 
the centres of galaxies to ones around protostars in star-forming 
regions, magnetic fields are almost always present and playing 
some role. With this paper, we begin a study of the properties 
of accretion discs around magnetised neutron stars. In particu- 
lar, we are interested in two kinds of system: X-ray pulsars and 
old neutron stars in the process of being spun-up (recycled) to 
become millisecond pulsars (MSPs). 



In X-ray pulsars the accreted matter transfers angular mo- 
mentum to or from the neutron star causing the spin frequency 
to increase or decrease at rates that are often hundreds of times 
faster than the typical spin-down rate in radio pulsars. Some of 
them are observed to be continuously speeding up or slowing 
down (with occasional reversals in these trends) while others 
show either little change in period or display erratic spin-down 
and spin-up behaviour (see, for example, the review by Bildsten 
et al. I1997I I. Exactly why the X-ray pulsars show such varied 
spin behaviour is still not clearly understood, but the magnetic 
field is almost certainly playing an important role in this. X-ray 
pulsars typically have magnetic fields of ~ 10'^ G and rotation 
periods in the range 10^ - 10"^ s. 

Old pulsars in the process of being recycled have lower mag- 
netic fields (typically ~ 10** G) and are being spun up to millisec- 
ond periods by means of accretion from their binary companion 
via an accretion disc. Their relatively weak magnetic fields al- 
low the inner edge of the disc to be close to the surface of the 
star. There has been recent evidence that the standard evolu- 
tionary model cannot explain the evolution of all MSPs, espe- 
cially the young ones with relatively high magnetic fields, e.g. 
PSR B1937-H21. Kiziltan & Thorsett (2009 ) showed that diff^er- 
ent MSPs must form by at least two distinct processes, but the 
nature of the second process remains unknown. 

Here we focus mostly on these "recycled pulsars", for which 
the distortions of the magnetic field are larger and occur at 
smaller distances from the central object, making the effects eas- 
ier to see. The behaviour for the X-ray pulsars is expected to be 
mainly similar, although one must then scale quantities because 
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of the larger inner radius of the disc and intensity of the magnetic 
field. 

One of the first important theoretical studies of accretion 
onto magnetised neutron stars was made by Ghosh, Lamb & 
Pethick ( 119771 1 who investigated the flow of accreting matter 
and the magnetic field configuration in the region inside the 
Alfven surface (taken to roughly coincide with the magneto- 
spheric boundary). In this pioneering work they were able to 
estimate the total torque exerted on the neutron star by the accre- 
tion disc, but at the cost of making many strong and questionable 
assumptions. For example, they assumed that the disc is com- 
pletely screened from the field except in a very small transition 
region, whereas there are some mechanisms (Kelvin-Helmholtz 
instability, turbulent diffusion and magnetic field reconnection) 
that allow the magnetic field to thread the disc across a very large 
region, as the same authors noted in a subsequent paper Ghosh 
& Lamb ( I1979al hereafter GL). Also, they used an ad hoc pre- 
scription for the azimuthal field, which was later shown to be in- 
consistent with having a stable disc beyond the corotation point 
(Wang 1987). They made a number of other assumptions which, 
while being approximations, are still used nowadays, such as 
having axisymmetry, a steady flow, infinite conductivity and a 
thin disc. The GL model was then used for calculating the torque 
in Ghosh & Lamb ( |1979bl l, where they pointed out that there are 
both positive and negative contributions to the torque, and that 
if the star rotates fast enough (with respect to the matter in the 
disc) the net torque can even be negative. 

In considering mechanisms for angular momentum loss from 
neutron stars, one should recall that there are at least two other 
processes which may well play a significant role: equatorial 
ejection of matter from the disc - this is the propeller mecha- 
nism (Davidson & Ostriker [19731 Illarionov & Sunvaev ll975l 
Shakura 1975); and gravitational wave emission, which would 
occur if the neutron star has some mass (or mass-current) asym- 
metry around its rotation axis (see Abbott et al. ( 12007l l). 

A significant step forward in studying these magnetic torques 
was made by Wang (fT987] l and Campbell ([19871 fT992l) . who re- 
considered the Bfj, generation mechanism. They succeeded in ob- 
taining an analytical expression for the azimuthal component of 
the magnetic field, which is still used despite the simplified na- 
ture of the model assumptions (e.g. having an exactly dipolar 
poloidal field and considering only the equatorial plane). In these 
models the induction equation is solved so as to find a steady so- 
lution for an axisymmetric magnetic field. They found that the 
generation of toroidal field is due only to the vertical gradient of 
the azimuthal velocity v^. If one further assumes that the disc is 
Keplerian (Q - Qk) and that the magnetosphere just above the 
disc is corotating with the star (with angular velocity Q.*), then 
one obtains that Z?^ oc Q^: - ^* ■ The magnetic torque is then cal- 
culated and, depending on the radii of the inner edge of the disc 
and of the corotation point r^, can either spin-up or spin-down 
the neutron star 

Our aim in the present work is to improve on this analysis 
and to find a stationary configuration for the magnetic field in- 
side the disc and in the suiTounding corona (the latter taken as 
being a layer above and below the disc), without making any 
leading order expansion or vertical integration, i.e. using a fully 
2D model. We want to study the effects of the relevant physical 
quantities (the radial, vertical and angular velocities, the diffusiv- 
ity, the accretion rate, etc.) on the structure of the magnetic field, 
and in turn on the magnetic torque, in order to give an estimate 
for the field structure in the disc which is more realistic than the 
analytic ones given by Wang ( 11987! ) and Campbell ( 11987! [T992b . 
We start by making use of the kinematic approximation (i.e. we 



consider the induction equation for a given profile of fluid quan- 
tities and do not solve the Euler equation to incorporate back 
reaction on the fluid from the field). Subsequently one can refine 
this model by including other important aspects, such as mag- 
netic back reaction. 

Miller & Stone d 19971 1 used a 2D model, including the ef- 
fects of non-zero resistivity in the MHD equations. They consid- 
ered the time evolution of an initial magnetic field profile but did 
not follow the evolution for long enough to reach a steady state. 
They found that regardless of the initial field structure (they used 
3 configurations) there was a rapid evolution of the disc, driven 
by angular momentum transport. In most cases, equatorial accre- 
tion results, either because accumulation of matter makes the gas 
pressure exceed the magnetic pressure or because the magnetic 
field geometry is such that polar accretion is inhibited. Polar ac- 
cretion only occurs when a strong global vertical magnetic field 
is included. Their simulations also confirmed the failure of total 
screening of the magnetic field from the accretion disc. 

Many investigators after Miller & Stone have simulated 
the magnetosphere-disc interaction, solving the full set of the 
MHD equations (see e.g. Romanova et al 2002 and KuUcarni & 
Romano va !2008! l . The work that we are caiTying out here should 
not be seen as being in competition with these analyses, but 
rather as being complementary to them. Our approach here is to 
use a succession of simplified models, becoming progressively 
more sophisticated, and to proceed step by step so as to fully un- 
derstand the effect of each successive additional feature as it is 
introduced. In large-scale numerical work, one sees the results of 
an interacting set of inputs within the scope of the adopted model 
assumptions and numerical techniques. Deconstructing this, so 
as to have a clear conceptual understanding of the role of each of 
the different components, remains a valuable thing to do and an 
approach which needs to be carried on alongside the large-scale 
simulations. The conceptual papers from the 1970s and 1980s, 
mentioned above, continue to be widely quoted and used as the 
basis for new research (see, for example, Kluzniak & Rappaport 
12007!) and our work here stands in the tradition of refining these 
approaches. 

Agapitou & Papaloizou (!2000!) looked for steady-state ax- 
isymmetric configurations of a force-free magnetosphere using 
simplified numerical calculations. However they focused on the 
global structure of the whole magnetosphere, rather on the de- 
tails of what happens inside the disc. Moreover in their model 
the disc is assumed to have only an azimuthal velocity and to 
be thin so that radial derivatives can be neglected with respect 
to vertical ones. They found that even under these circumstances 
the poloidal field can differ significantly from the dipole field 
of the central star, and that the magnetic torque can be much 
smaller than that estimated assuming B, ~ Bdip- 

Following a somewhat similar strategy, we aim to make im- 
proved calculations by using a more general profile of the veloc- 
ity field (with all of the components being allowed to be non- 
zero); by not neglecting any of the derivatives present in the in- 
duction equation; and by using a better resolution (i.e. a finer 
grid) and a much larger numerical domain in the radial direc- 
tion so as not to have the results being dependent on the outer 
boundary conditions. 

The plan of this paper is as follows. After the present 
Introduction, the details of the model adopted here are given in 
Section |2] while the equations solved are presented in Section 
[3] As we will show there, since the system is taken to be ax- 
isymmetric and in a stationary state, the induction equation can 
be split into two parts, one in which only the poloidal compo- 
nent of the magnetic field appears and the other containing all of 



L. Naso and J.C. Miller: An investigation of magnetic field distortions in accretion discs around neutron stars 



3 



the components. In this paper we present results for the poloidal 
field; the toroidal one will be discussed in a subsequent paper 
The numerical code used is described briefly in Section [3] and 
then in more detail in the Appendix, where we describe both 
the algorithm used and the tests which we have performed for 
verifying it. In Section |4] we present our results and Section |5] 
contains conclusions. 



2. Our Model 

In this study we consider disc accretion by a neutron star having 
a dipolar magnetic field with magnetic moment ^. We assume 
that the star is rotating around its magnetic axis, and that this axis 
is perpendicular to the disc plane; also, we assume that the fluid 
flow is steady and has axial symmetry everywhere. Of course, all 
pulsars must be oblique rotators in order to be observed as such, 
and so the system would not then be axisymmetric as a matter 
of principle. However we are not concerned here with studying 
the aspects of accreting systems affected by the inclination angle 
(e.g. the impact region or stellar precession) and the principles of 
the interaction between the magnetic field and the plasma in the 
disc are not affected by the inclination. Therefore we feel that re- 
taining the assumption of the magnetic axis being perpendicular 
to the disc is satisfactory for our purposes. For our calculations, 
we use spherical coordinates (r, 9, 4>), with the origin being at 
the centre of the neutron star 

The general flavour of our model is similar to the GL one. 
We suppose that at large radii the magnetic pressure is negligi- 
ble with respect to the gas pressure and the disc can be described 
as a standard o'-disc. As one moves inwards however the mag- 
netic field becomes progressively stronger, and may eventually 
dominate over the gas pressure so that most of the matter leaves 
the disc following magnetic field lines and the disc is disrupted. 
We note that some equatorial accretion continues to be possible 
even with a rather high stellar magnetic field, as shown by Miller 
& Stone flWn . 

There are some differences between our model and that of 
GL. We are not assuming the existence of a wide transition zone 
where the field is progressively screened, eventually disappear- 
ing at the outer boundary of the zone. Also, we do not force the 
field to be zero at the outer boundary; we instead take it to be 
dipolar there and we do the same at all of the other boundaries 
as well. We put the radial outer boundary very far away so that 
the solutions in the zone of interest, which goes from the inner 
edge of the disc out to the light cylinder, will not be very sen- 
sitive to the conditions at the outer edge. We note here that, if 
we take the innermost part of the disc to have very high diffusiv- 
ity, we naturally find a narrow inner region followed by a broad 
outer one, however the behaviour of the field in these regions is 
quite different from that in the GL model, as will be shown in 
the results section. 

As regards the inner edge of the disc, its precise location is 
still an open issue and several prescriptions have been suggested, 
however none of them differs very much from the Alfven radius 
calculated using a dipolar magnetic field (i.e. the radius where 
the pressure of the dipole field equals the gas pressure). For our 
present model we take this as giving the inner edge of the disc; 
for subsequent models, a possible improvement will be to cal- 
culate the Alfven radius using the stationary magnetic field ob- 
tained in this work instead of the dipolar one. 

We suppose that all around the pulsar there is vacuum, ex- 
cept for where we have the disc and the corona, and that the 
field remains dipolar from the surface of the star until it reaches 



the matter in the corona. In reality, taking vacuum is not com- 
pletely correct, both because the density in the magnetosphere is 
not zero and because between the star and the disc there is the 
matter which is accreting onto the neutron star. For the latter, we 
suppose that the flow of this material is perfectly aligned with 
field lines and that it has very low density, so that we can neglect 
its influence on the magnetic field structure. Furthermore we are 
introducing a low density corona in order to have a transition 
zone between the disc and the vacuum. We also allow the veloc- 
ity and the diffusivity to have different values in the corona from 
those in the disc (we will comment on this later in the present 
Section). 

Summarizing, we model the system as being composed of 4 
regions (see figure [T]i: (1) a central object, surrounded by (2) an 
accretion disc, on top of which there is (3) a corona; all of the 
rest is taken to be (4) vacuum. Each physical quantity is allowed 
to have different values in each of the regions. Our numerical 
domain covers regions (2) and (3). Since these two regions are 
surrounded by vacuum and the stellar field is dipolar, we im- 
pose dipole boundary conditions at all of the boundaries. We are 
aware that in real astrophysical systems this dipole condition at 
the boundaries can be rather drastic since the region outside the 
neutron star is not vacuum and the field is distorted well before 
reaching the disc. However at our level of analysis the results are 
not very sensitively dependent on the precise profile chosen for 
the magnetic field at the boundaries, and we are here focusing on 
studying the influence on the magnetic configuration of the ve- 
locity field and diffusivity. That this is reasonable is confirmed 
by the fact that the magnetic configuration which we obtain is 
rather similar to that given by the numerical simulations of other 
authors (Miller & Stone, 1997), which had a different treatment 
of the boundary conditions. We will comment more on this in 
Section m 

We take the ratio H/r to be constant (where H is the semi- 
thickness of the disc), and so the entire upper surface of the disc 
is located at a single value of the colatitude (as also is the case 
for the corona). In particular, we take an opening angle of 8° 
for the disc (measuring from the equatorial plane to the top of 
the disc), implying H/r ~ 0.14, and of 10° for the disc plus the 
corona. 

As will be described in Section 3, in this paper we consider 
only the poloidal component of the induction equation. Once the 
magnetic diffusivity rj and the poloidal velocity Vp are specified, 
this can be solved to obtain the configuration of the poloidal field 
without entering into details of the toroidal component. As re- 
gards the turbulent magnetic diffusivity, we take this to have a 
constant value t/o in the main disc region and to be 77c = 100 ■ 770 
in the corona and in the inner part of the disc (we join the differ- 
ent parts smoothly, using eiTor functions). We take higher values 
in the latter regions because the density is lower there and so we 
expect the effects of turbulence to be enhanced. However as we 
move away from the corona into the vacuum the density drops 
to zero and turbulence eventually disappear^]. 

As regards the velocity: for vv we use the expression given 
for the middle region of ff-discs^lby Shakura & Sunyaev (1973); 
while we put vg to zero inside the disc, although allowing it to 
be non-zero in the corona as necessary in order to be consistent 
with the dipole boundary conditions. In Section l3T2l we find that 



' It is only the turbulent diffusivity rjj which disappears, the Ohmic 
one ?7ohm always be present, therefore in vacuum 77 = Tyohm ^ '7T- 

' For the parameter values which we are using, both the inner edge 
of the disc and the outer boundary of our domain of interest lie within 
the "middle region" of the a-disc model. 
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v^Br - VyB^ + -drirB^) 



(5) 



r sin0 



dg{B^ sin 6) 



The first two of these equations have the same expression inside 
the large square brackets and we call this T. We can then reduce 
these two equations to a single equation and write the entire sys- 
tem as: 



Fig. 1. Geometrical representation of our model (the drawing is 
not to scale). The radius values which we use were: ~ 10 r^, 
R,r ~ 22 Tg and Ri^ ~ 1 15 /"g and the opening angles are 8° for 
the disc alone and 10° for disc plus corona. The numerical outer 
boundary is much further out than the region of interest shown 
here; the grid continues until r ~ 750 (this is the reason for 
the dashed lines). 



whenever a dipole field is a solution of the induction equation, 
a precise relationship must hold between v, and vg, subject to 
certain reasonable conditions. We use this relation to calculate 
Vg in the corona. 

Figure[T]shows a schematic representation of the geometry of 
our model: the corona and the two parts of the disc being shown 
(the transition region is thin and around r = R,r). The dashed 
lines indicate that in solving the induction equation we use a 
much larger numerical domain, so that in the zone of interest 
the solution is almost independent of the location of the outer 
boundary. 



3. The Equations 

The induction equation relates the time variation of the magnetic 
field to the properties of the plasma. If one considers mean fields, 
it can be written as: 



r sin 6 



5,B = V X (v X B + £ - 77ohmV x B) 



(1) 



where ?7ohm - c^/^-ncr is the Ohmic difFusivity and £ is the 
turbulent electromotive force, £ = {\' xB'), with the primes 
denoting small-scale turbulent quantities. A common procedure 
is to expand £ in terms of the mean field and its derivatives 
and within the first order smoothing approximation one has 
£ = aB - 77tV X B, where the o-B term generates the so-called 
Q'-eff'ect that is fundamental for having a dynamo action. In the 
present paper we neglect this effect, which could however be in- 
cluded in a more elaborate subsequent model. This is in line with 
our approach, which is aiming at understanding, one at a time, 
the effects of the various elements which characterize the system 
of an accretion disc around a magnetised neutron star. 
The induction equation then reduces to: 



(9,B = V X (v X B - 77V X B) 



(2) 



where - T/ohm + Jfl ~ because the turbulent diffusivity is 
much larger than the Ohmic one in the disc and in the corona. 

Since we are interested in stationary solutions, we put the 
time derivatives to zero, as we do also for the <p derivatives since 
we are assuming axisymmetry. Then in spherical coordinates the 
three components of equation ^ are: 



- dglsint 



^ dr^r 



VrBg - VgBr [drirBg) - dgBr] 

r 

VrBg - VgBr [dr{rBg) - dgBy] 

r 



(3) 
(4) 



Q ^ dr\r 



pBr - VrB^ + -drirBg) 



(6) 
(7) 



dg I vgB^ - v^Bg l—:dg{B^ sin 0) 



where is a generic constant. 

We notice that T depends only on the turbulent magnetic dif- 
fusivity rj and on the poloidal components of the velocity field 
and r and component of the magnetic field. It is clear therefore 
that equation (|6]l alone governs the poloidal part of the magnetic 
field and is independent of any azimuthal quantity, while to ob- 
tain the toroidal component of the magnetic field one has to solve 
the further equation 

In this paper we concentrate only on solving for the poloidal 
field, using equation (|6]l, for which no knowledge of behaviour 
in the direction is required. In a forthcoming paper we will use 
the results obtained here to solve equation (|7]) and calculate the 
toroidal field component. 

In order to guarantee that the condition V ■ B = is satisfied 
and to be able to calculate the magnetic field lines easily, we 
write the magnetic field components in terms of the magnetic 
stream function S, which is implicitly defined by the following 
two equations: 



Br = 



Bg = - 



1 



• sin 6* 
1 

rsmO 



dgS{r,ff) 



drSir, 0) 



(8) 
(9) 



With this definition, the axisymmetric field B is always 
solenoidal and isolines of S represent magnetic field lines. 
Substituting these expressions into equation (|6]l, we obtain an 
elliptic partial differential equation (PDE) for S in terms of r 
and 0: 



o It / COt0 Vfl \ Vr k 

d].s + ^dls-[^ + —\dgS--drS^- (10) 

\ rrjj rj rj 

where k is the constant introduced in equation (|6]) and Vy, vg and rj 
are non-constant coefficients. This equation can be solved once 
boundary conditions and values for the coefficients have been 
specified. 

3. 1 . Velocity field and turbulent diffusivity 

In Section [2] we qualitatively described the profiles of velocity 
and diffiisivity that we are using. Here we give the precise ex- 
pressions. 

For the velocity, we use the expression given by Shakura & 
Sunyaev (fT973T l: 

Vr(r) = 2 lO^-a^'^-in^^^-m-^'^i3/rf'^ [l - (3/r)°-^]"^^%ms"'(ll) 

where the radius r is expressed in units of the Schwarzschild 
radius, m is given in solar mass units, m is in units of the critical 
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Eddington rate and a is the standard Shakura-Sunyaev viscosity 
coefficient. Using typical values (a = 0.1, ih = 0.03 and m - 
1.4) one obtains: 



Vr(r) « 7.3 10"* ■ [l - ^^^ms"' 



(12) 



For the other component of the poloidal velocity, vg, we set 
this to zero in the disc and use a non-zero profile in the corona. 
The formula for this is calculated in Subsection l3.2l we antici- 
pate here the result: 



Vg 



in the disc 

1 V, tan 6 in the corona 



(13) 



For the diffusivity, we first construct two auxiliary functions, 
rjg and rjr, giving the profiles along the 6 and r directions respec- 
tively: 



7]r = 



1 -i-erf 
1 + erf 



-r + r„ 



d. 



1 + erf 



A }] out 



(14) 
(15) 



where 6c is the colatitude of the upper surface of the disc, r^^in 
is the radius of the boundary between the inner region and the 
main disc regiorQ and dg and dr are the widths of the error func- 
tions used for the angular and radial profiles respectively. We 
then combine equations (fT4l i and ( fTsT i to get the global rj: 



i](r, 6) = 7/0 



i + + ri.) ■ (- - 1 
2 \m 



(16) 



where 770 and ?7c are dimensional quantities (with units of cm 
s"'), the former giving the value of the diffusivity in the main 
disc region and the latter giving the value in the corona and the 
inner disc region. 



3.2. Dipolar solution, an analytic constraint 

In order to obtain a profile for vg we consider the situation when, 
from the top surface of the corona (at 6 = 0surf) down to some 
6-6, the stationary magnetic field is dipolar, i.e. (Br, Bg, B^) - 
^2/i_rasfl^ /ijmfl^ where ju is the magnetic dipole moment. Since 
a magnetic dipole field is current-free (i.e. J = V x B = 0) 
equation ^ becomes: 



V X (v X B) = 



(17) 



Following a procedure similar to that used for obtaining equa- 
tions (|6]l and Q, we go to spherical coordinates, write out the 
three component equations and group the poloidal terms. This 
gives: 



(18) 
(19) 



V,. Bg -VgBr - : — - 

r sin0 

dr{r V0 Br) = -dg{v^Bg) 



where is a generic constant. In order to calculate it we consider 
a path with constant 6, e.g. 9 - 6* e [flsurf,^]; along this path 
equation (fTSl ) gives: 



/u sins'* \ ^ J I k 

Vr I — \- Vg {ji cos 6 ) - r 



sinfl* 



(20) 



^ r^^ou, appearing in the expression for 77,- is located shortly before To, 
and far away from the zone of interest. 



where all of the terms in the brackets are constant. To investi- 
gate further the properties of this equation, consider the situa- 
tion if Vg - 0. Equation (l20l i would then imply that either Vr oc p- 
(which is not reasonable for an accretion flow) or = so that 
there is no accretion at all. This implies that if we have an accret- 
ing flow in a region with a dipole magnetic field then vg must be 
non zero. We do not know the exact profiles of v,- and vg, but it 
is not plausible that the left hand side increases as p- and so we 
need to choose ^ = 0. Therefore from the last equation we get: 



sin0 



Vr ■ 



— Vg COS f 



1 

Vg — — tanflvr 



(21) 



Equation (1211 1 implies not only that if there is a non-zero radial 
velocity then there must be a non-zero vertical velocity as well, 
but also that the vertical velocity is larger than the radial one (for 
= 81° one has ~ 13 Vr). 

The behaviour of the azimuthal velocity v^ can be investi- 
gated using equation (fT9] l. As expected = is a possible so- 
lution but V0 = VKep is not, meaning that having a dipolar field is 
not consistent with having a Keplerian angular velocity profile, 
whereas it is consistent with having no rotation. It is interesting 
to notice that there are also some non trivial profiles which are 
solutions. For example v^ oc r^^ ■ sin^^^"^ 9, where 5 is a posi- 
tive integer, gives a set of possible solution^ Therefore if one 
supposes that v^ decreases with r as a power law, then it must 
have also a dependence on 9 in order for the magnetic field to 
be dipolar We recall that in the works of Campbell and Wang it 
is precisely the vertical gradient of the angular velocity that pro- 
duces the toroidal field. Here we have shown that it is possible 
to have a non-zero vertical gradient of the angular velocity and 
still have a zero toroidal magnetic field. However, we stress that 
we are not giving physical explanations for having these kinds 
of velocity profiles. 

We note that equation ( fTTl i has been solved in the context of 
stellar winds by Mestel (1961); his result was that the poloidal 
magnetic field and velocity field need to be parallel. Our result 
that VrlBr - vg/Bg (from equation (fTST l with = 0) is in agree- 
ment with this. 

3.3. Solution metliod 

Before solving equation ( fTOl i we put it into a dimensionless form, 
by scaling the quantities in the following way: 

r = rrg 6^6 S = SSo (22) 

where the hat quantities are dimensionless, is the 
Schwarzschild radius, 770 is the value of the diffusivity in the 
main disc region and So is a reference value for the stream func- 
tion, calculated as the value for a dipolar field on the equator of 
the neutron star Substituting into equation (fTOl i we get: 



'g 'g' 



'cote 

2 '7 



cot 9 Vg Tg 



Vr So 

SodgS d-rS 



Vrr. 



dnS -d-rS = 



0(23) 



(24) 



where we go from (1231 ) to (l24l i by dividing both sides by So/r^. 
We rename the variables (r = x and 9 - y) and obtain the fol- 
lowing dimensionless form for the equation: 

T - It' / cot V VgrA , Vrr„ 

dls + —dlS - — f + — \d,S- — d,S = 

- \ XT] I ' T] 



(25) 



This can be seen by writing as the product of a function of r and 
a function of 6 and carrying out some simple algebra. 
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where [v^] - [ve] - cm s"', [rg] = cm, [rj] 
Vr rg/rj and vg rg/77 are dimensionless coefficients. 

As mentioned above, the reference value for the stream func- 
tion is taken to be its value at the equator of the neutron star 
If we call the star's radius ro rg then So = S'^^^iro r„, njl) - 
/^('^Cg) The solution of equation (IZST i gives the dimension- 
less stream function S. In order to calculate the magnetic field, 
we first calculate B^o\ using the dimensionless versions of equa- 
tions (O and (|9]l, and then multiply by Bq to go to physical units, 
where Bo = B'^g'^{r - ro rg, 9 = n/2). 

We have solved equation dZSl l with the Gauss-Seidel relax- 
ation method, which uses a finite difference technique, approx- 
imating the operators by discretizing the functions over a grid. 
The value of the stream function at any given iteration step is 
written as a function of its value at the previous step, or at the 
same step but at a location which has already been calculated. 
Details of the numerical scheme are given in the Appendix. 

Before applying the numerical scheme to the real problem 
that we want to solve, we performed a series of tests on the code, 
which are described in detail in the Appendix. We used sev- 
eral configurations, with many different numbers of grid points, 
profiles for the velocity and diffusivity, initial estimates for the 
stream function, locations for the outer radial boundary of the 
grid and values for the iteration time step. In this way we have 
checked the stability and convergence, have optimized the iter- 
ation procedure and have determined where to place the outer 
radial boundary of the grid (which needs to be far enough away 
so that the outer boundary conditions do not significantly influ- 
ence the solution in our region of interest). 

All of the results presented in the next section have been ob- 
tained using a 1000 x 20 grid and with 2^'^ ~ 1.7-10^ total 
iterations. With these settings, we have always obtained residu- 
als of the order of lO"''* or less, and the resolution is Ax ~ 0.7 rg 
and A0 -0.5°. 



4. Results 

From equation (IZST l. we see that in order to calculate the stream 
function S, and hence the poloidal magnetic field, we have to 
specify vv(r, 6), vg{r, 6) and ri(r, 6). For the vertical velocity we 
follow the prescription given by equation (1211 1: the profiles used 
for the other two functions have been described previously (see 
equation (fTTb for Vr and Section |2] for rj). Here we describe how 
the magnetic field configuration changes when we modify these 
two functions. 

In figures |2]and[3] we show the magnetic field lines calculated 
with four values of vo, i.e. with different accretion rate^ For fa- 
cilitating the comparison we show also a dipolar magnetic field 
(dotted). The field lines are labelled with the radial coordinate 
where the dipole field imposed at the top boundary would cross 
the equatorial plane if not distorted. We can see that if vv were 
zero, the field would not be distorted at all from the dipolar con- 
figuration and increasing the velocity then creates progressively 
more distortion. The degree of distortion depends on the location 
in the disc: in the inner part, where the field is strongest, it is most 
able to resist distortion; further out, the field is weaker and it be- 
comes progressively more distorted. From the top panel of figure 
[3] we can see that the distortion increases with distance. For ex- 
ample if one considers only the configuration with vo - 10"^, then 
one sees that the distortion on the equatorial plane at r = 15 rg 



Vo = 0, 77o = 10'", R„ = 



' For the Shakura-Sunyaev model, the radial velocity is proportional 
to m-l^ if the mass of the central object and the viscosity a are kept 
fixed. 




= 10 , 7]o = 10 



Numerical solution 
Dipole 



20 



30 



40 



50 



Fig. 2. Magnetic field lines from the numerical solution (solid) 
and those for a dipole (dotted) for comparison. The top panel 
refers to a case with no accretion, where the field remains exactly 
dipolar, while in the bottom panel vq = 10"* cm s and the field 
is distorted. The diffusivity t/q has the same value in both panels 



(10' 



'). If we use 770 = 10" cm^ s ' and vo - 10^ cm s 



we obtain the same result as shown in the bottom panel (/?„ 
for both). 



is roughly zero, at r = 25 rg it is slightly larger than 1 rg and at 
r = 45 rg it is about 5 rg. 

According to the behaviour of the magnetic field lines, we 
can divide the disc into three regions: (1) an inner region, where 
the lines are not distorted very much away from the dipole; (2) a 
main region, where the distortion is very large and (3) the region 
in-between the two, which we call a transition region, where 
there is an accumulation of field lines. This means that in the 
transition region there is an amplification of the magnetic field 
(see figure|4|i. 

In addition to varying the radial velocity, we also consider 
the role of the diffusivity. The results show that when we change 
770, we get the opposite behaviour to that seen when varying the 
velocity, i.e. a larger tjq gives a smaller distortion. Actually what 
really matters is not the velocity or the diff'usivity alone but their 
ratio. This is not surprising since in the equation which we are 
solving (equation ( l25T l) the quantities only appear in this ratio 
(bearing in mind that vg is taken to be either zero or propor- 
tional to V,.). In fact the magnetic Reynolds number 7?,^, which 
describes the general solution of the induction equation Q, is 
built from them: it is defined as = Iq ■ vq/jjo, where Iq, vq 
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Fig. 3. The same as in figure |2] but with different values of vq. 
The top panel is for a case with vq = 10^ cm s"' while the bottom 
one is for vq - 10'' cm s"'. If we use t]q = 10^ cm^ s"' and 
vo - 10^ cm s"', we obtain the same as in the bottom panel 
(Rn, = 400 for both). 



and ?7o are respectively a characteristic length, velocity and dif- 
fusivity. This parameter gives the relative importance of the two 
terms on the right-hand side of the induction equation. For large 
Rm, we are in the regime of ideal MHD with the magnetic field 
and plasma being frozen together, for low instead, the field 
and plasma are almost decoupled and the field simply diffuses. 
In accretion discs the radial velocity is usually many orders of 
magnitude smaller than the azimuthal velocity. In our case the 
Reynolds numbers calculated using the two velocities differ by 
about five orders of magnitude, if one takes the Keplerian veloc- 
ity as the characteristic velocity. For our present calculations 
only the radial motion is relevant, since = in the disc and 
is proportional to v, in the corona, while does not appear in 
the equation that we are solving now. The value of R^n reported 
in figures|2]and[3]is therefore the one calculated taking the char- 
acteristic velocity to be the radial velocity. This is a key point 
for the present considerations: the relevant magnetic Reynolds 
number is that calculated with the radial velocity and not that cal- 
culated with the azimuthal velocity (which gives a much larger 
value). 

The panels in these figures are clearly showing that the dis- 
tortion of the field is proportional to the magnetic Reynolds num- 
ber calculated in this way. This happens because with increasing 
Rm the freezing condition gets progressively stronger so that any 



fluid motion perpendicular to the magnetic field lines encounters 
more and more resistance. Therefore, since the velocity field is 
fixed, the magnetic field has to change. Figures |2] and [3] not only 
show that modifications in the magnetic field lines increase with 
Rjn, but also that their shape is consistent with what is expected 
from considering the flux freezing condition in the case of a con- 
ical flow (which is what we have in the disc). 

However the actual value of the magnetic Reynolds number 
is somewhat arbitrary, because in general there is no unambigu- 
ous way of defining the characteristic length, velocity and diffu- 
sivity of a given system. In our case we choose /q to be the radius 
of the inner edge of the disc, vo to be the radial velocity at the 
inner edge of the disc and 770 to be the value of the diffusivity in 
the main disc region. One could also make a different choice for 
the characteristic length /q, such as taking this to be the radius 
of the star or the average height of the disc; the trend of having 
larger distortions for larger values of would of course be seen 
in all cases, but the switching on of the distortions would occur 
at different threshold values of 7?^. 

We have already noted that the distortion varies with posi- 
tion, and so it is clear that a single global parameter cannot give 
a sufficiently detailed description in all parts of the system. It 
is therefore convenient to introduce a new quantity which we 
call the "magnetic distortion function" We define this in the 
same way as the magnetic Reynolds number but, instead of tak- 
ing characteristic values for the velocity and diffusivity, we take 
the local values: 

Dm(r, 0) = , „ (26) 

T]{r, 6) 

This function gives the relative importance of the two terms on 
the right-hand side of the induction equation at every point of the 
disc, rather than giving just a global measure as with the standard 
magnetic Reynolds number We then expect the advection term 
(V X (v X B)) to dominate in the regions where > 1 and the 
diffusion term (V x (77V x B)) to dominate when < 1 . This 
then explains why we find three regions inside the disc: the inner 
region corresponds to the zone where D,n <g: 1, the main region 
to Dm » 1, and the transition region to intermediate values of 
Dm. This correspondence is made clear in figure 21 where we 
show the component of the magnetic field, the dipolar profile 
and the magnetic distortion function. We recall that the jump in 
Dm follows from the profile chosen for 77, i.e. we use a larger 
value of the diffusivity in the inner part of the disc. 

Another important aspect of the magnetic distortion function 
is that the degree of arbitrariness in its definition is smaller than 
for the standard magnetic Reynolds number, since it is defined 
using only one characteristic value, Iq. In addition there is a quite 
natural way for choosing Iq. By looking at equation dZSl ) one can 
see that, if Iq = rg, the magnetic distortion function is akeady 
there in the equation (it is the coefficient of the partial derivative 
of S with respect to x). The value rg is coming from the way in 
which we are scaling the lengths. If we had chosen a different 
unit for the lengths, say T, then we would have had to choose 
Zo = r if we wanted Dm to appear directly in equation dZSb . We 
can then think of /o as a quantity needed to make the ratio v/rj 
dimensionless, and the most natural choice for this is the char- 
acteristic scale being used as the unit length. 

Summarizing, we can describe the magnetic field configura- 
tion in the accretion disc by saying that magnetic field lines that 
enter the disc in the main region (Dm » 1 ) are pushed towards 
the central object, whereas those which enter the disc in the in- 
ner region (Dm 1) are almost unmodified. The result is that 
in between these two regions there is an accumulation of field 
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Fig. 4. The magnetic field Be as obtained from the numerical 
simulations for the configuration with - 40 (dotted) com- 
pared with those for a perfect dipolar field (dashed) at the equa- 
torial plane. The solid line is the magnetic distortion function. 
The scale for is shown on the right vertical axis. The mag- 
netic field is measured in units of the stellar field, Bq = 3 10** 
G, as shown on the left vertical axis. 



lines, and so there is an amplification of the magnetic field there, 
as can be seen in figure |4] 

In order to test the sensitivity of the results to the boundary 
conditions, we have experimented with several different profiles 
of the magnetic stream function outside the disc. We have used 
three additional profiles: one which gives a magnetic field with 
spherical field lines, one which gives a magnetic field with ver- 
tical field lines and another one which gives field lines inclined 
at an arbitrary angle to the vertical axis. We have chosen these 
magnetic field profiles by comparison with the results from the 
simulations by Miller & Stone (I1997I I. For magnetic field lines 
entering the disc at the same locations, their shape within the 
disc varies hardly at all in the different cases, although the field 
intensity does vary significantly. However, we are focusing here 
on the shape, looking at the distortion of the field lines, which 
does not depend sensitively on the boundary conditions used and 
is instead mainly governed by the magnetic distortion function 

ZJm. 

In order to understand better the influence of the magnetic 
distortion function on the magnetic field structure, we varied 
and saw how the field changed. We used three new profiles for 
Dn, different from the previous one which we then considered 
as a reference. In the first profile we increased the value of 
in the inner part of the disc and left the rest unmodified, in the 
second one instead we lowered in the main region of the 
disc and did not change the inner part, and in the last one we just 
changed the width of the transition between the low and high 
values of D^. We then calculated the poloidal magnetic field 
and the results are presented in figure |5] where the top panel 
shows the different profiles of the magnetic distortion function 
and the bottom one shows the component of the magnetic field, 
referring to the equatorial plane in both cases. 

Considering this figure, we can summarize the influence of 
the magnetic distortion function with four comments: (1) chang- 
ing the value of D,,, in the small inner part by a factor of 5 leaves 
the magnetic field almost unchanged, (2) on the other hand, the 
magnetic field is very sensitive to the width of the transition in 
Dm and to its value in the main region, in particular (3) the po- 
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Fig. 5. Top panel: magnetic distortion function in the equatorial 
plane. Bottom panel: 9 component of the magnetic field in the 
equatorial plane. For both panels: the solid line shows results 
from the previous analysis (in the bottom panel the dipolar field 
is shown with a thick solid line); the dotted lines refer to profile 
1 (larger value of in the inner region only); the dashed lines 
refer to profile 2 (lower value of in the main region only); the 
dot-dashed lines refer to profile 3 (same values of D,^ in inner 
and main region, but the transition region is wider). 

sition of the peak in Bg is related to the width of the transition 
and (4) the deviations away from the dipole field are mainly gov- 
erned by the value of Dn, in the main region. We can go further 
and consider the radial derivative of D„^, which is shown in fig- 
ure |6] for all of the profiles used. From this we can see that the 
position of the peak of Bg is strongly connected with the posi- 
tion of the minimum in drD^, and that the maximum amount 
of magnetic distortion is related to the depth of the dip in the 
derivative of D^. 

Finally we comment on the behaviour of the magnetic stream 
function in the equatorial plane. In order for the magnetic field 
to have a local minimum or maximum, its first radial derivative 
must, of course, be zero. This condition can be written in terms 
of S as: 

d^,S - -drS = (27) 
r 

In the numerical simulations S is always decreasing with r and 
so the only way to get a local minimum or maximum is for S to 
go through a region where its concavity is reversed and becomes 
negative. An example of this is shown in figure|7] where we plot 
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Fig. 6. Radial derivative of the magnetic distortion function in 
the equatorial plane. The profiles refer to the same case as in 
figure |5] Comparing with the bottom panel of figure |5] one can 
see that the peaks in the magnetic field Be occur at almost the 
same locations as where dr £>m has a minimum, and that the am- 
plitudes of the distortions are proportional to the absolute values 
of the minimum of d,- D^. 
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Fig. 7. The magnetic stream function S for the reference profile. 
The dotted lines are drawn where S changes concavity, while 
the open circles mark the locations where Bg has a minimum or 
a maximum, i.e. r ~ 18 rg and r ~ 22 r^; they are in the region 
where the concavity of the magnetic stream function is negative. 



the stream function for the reference case. We can see that for 
r < 16 rg and r > 22.5 S has positive concavity, while be- 
tween these two values it is negative. The dotted lines in the fig- 
ure delimit the regions of positive and negative concavity, while 
the open circles are drawn at the locations where Bg has a mini- 
mum or a maximum. As expected these points are in the region 
of negative concavity. 



5. Conclusions 

In this paper we have begun a systematic study of the mag- 
netic field configuration inside accretion discs around magne- 
tised neutron stars, which is intended as being complementary 
to the large numerical simulations being carried out elsewhere. 
We have assumed that the star itself has a dipolar magnetic field. 



that is rotating around its magnetic axis and that this axis is per- 
pendicular to the disc plane. We have also assumed that the flow 
is steady and has axial symmetry everywhere. Our strategy was 
to start the analysis with a very simple model, where we made 
the kinematic approximation and solved the induction equation 
numerically in full 2D, without making any leading order ex- 
pansion. This initial model will subsequently be improved by 
including the magnetic backreaction on the fluid flow. 

We have shown that it is possible to separate the calculation 
for the configuration of the poloidal magnetic field from any az- 
imuthal quantities. This is a key point and has the consequence 
that the effective magnetic Reynolds number is that calculated 
with the radial velocity rather than the azimuthal one (which is 
very much larger). We have here considered only the poloidal 
component; the toroidal one will be addressed in a subsequent 
paper 

We have modelled the system as being composed of four re- 
gions (see figure[TJ: the central neutron star, the disc, the corona 
(taken to be a layer above and below the disc) and all of the 
rest, which is taken to be vacuum. We suppose that the stellar 
magnetic field remains dipolar until it reaches the corona. At 
that point it begins to feel the presence of the fluid flow and the 
magnetic field lines are pushed inwards, thus creating distortions 
away from the purely dipolar field. 

We have studied the response of the magnetic field to 
changes in the velocity and the diffusivity, finding distortions 
away from dipolar increasing with the radial infall velocity 
and decreasing with increasing diffusivity. The underlying be- 
haviour is that the distortions increase together with the magnetic 
Reynolds number R^n (which governs the flux freezing condi- 
tion) where the ratio vo/rjo appears. 

However a single value of cannot take into account any 
large changes in the magnitudes of the velocity and the diffu- 
sivity through the disc, since it is defined using single charac- 
teristic values. Therefore in order to have a sufficiently detailed 
description of the system, we have introduced a magnetic distor- 
tion function Dm = '"^"^^y based on local values of the quanti- 
ties concerned, so that in the regions where s> 1 or 1 
one should expect to have large or small distortions respectively. 
We expect the turbulence to be enhanced in the regions of lower 
density (the corona and the inner part of the disc), therefore in 
our model we use a larger value of 77 in these regions, giving a 
smaller value for D^. As clearly shown in the panels of figures 
|2]and[3] the disc can be divided into three parts: (1) an inner 
region, where 1 and the distortions are negligible; (2) a 

transition region, where is rapidly increasing and magnetic 
field lines accumulate; (3) a main region, where » 1 and the 
the distortions are very large. Considering D„i can be very useful 
when analysing results of large numerical simulations. 

Comparing our results with previous literature, we can con- 
firm the idea of dividing the disc into two principle regions: an 
inner part, where the magnetic field is strongest, and an outer 
part, where the magnetic field is weaker and gently decaying. 
However the behaviour that we find for the field in these regions 
is very different from that of the Ghosh & Lamb model ( 1979al) 
and we find it convenient to include a third zone, to be consid- 
ered as a transition between the two principle ones (see figure 
|4|i. In the inner boundary layer of the GL model, the magnetic 
field is reduced by screening currents by a factor of 80 per cent, 
while in our case the field is barely modified in the first region. In 
the transition region instead we see a completely new effect: the 
field is amplified and has a local maximum. The properties of the 
maximum (i.e. its location, the peak magnitude of the field and 
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its behaviour in the neighbourhood) are well described by the 
magnetic distortion function D^, in particular by its maximum 
value, by the width of the transition between the low and high 
value and by the behaviour of the radial derivative drD^ (i.e. by 
the location and the magnitude of its minimum). The behaviour 
in the outer zone is instead quite similar to that in the GL model, 
with the field decaying and being smaller than the dipole one at 
any give location. 

As regards the magnetic field geometry, our results resemble 
rather closely those obtained by Miller & Stone (il997^ (compare 
our figure[3]with the top panels of their figure 3), despite the fact 
that they solved the full set of MHD equations whereas we have 
solved just the induction equation and with different conditions 
at the top of the disc. Moreover we have found that the distortion 
of the field lines inside the disc depends very little on the pro- 
files outside it. We conclude that the distortion pattern seen for 
the poloidal component of the magnetic field should be rather 
a general result, related only to the fundamental aspects of the 
system, and that the distortion is not at all negligible for typical 
values of the accretion rate (see figures |2] and |3]l. 

The immediate next step will be that of calculating the 
toroidal magnetic field component, by solving equation (|7]i, and 
consequently the magnetic torque. After that we will move to 
the next model, by removing the kinematic approximation and 
solving the Euler equation as well as the induction equation. 
By comparing the results of that calculation with those for the 
present kinematic model, we will be able to understand the ef- 
fects of the magnetic backreaction, which we expect to affect 
mainly the toroidal magnetic field component. 

Acknowledgments 

It is a pleasure to thank Alfio Bonanno, Claudio Cremaschini 
and Kostas Glampedakis for stimulating discussions; this work 
was partly supported by CompStar, a Research Networking 
Programme of the European Science Foundation. 



References 

Abbott B. et al., 2007, Phys.Rev.D, 76, 042001 

Agapitou V. & Papaloizou J.C.B. 2000, MNRAS, 317, 273 

Balbus S.A. & Hawley J.F. 1991. ApJ, 376, 214 

Bildsten L., Chakrabarty D., Chiu J. et al. 1997, ApJS, 1 13, 367 

Bonanno A. & Urpin V. 2006, Phys. Rev. E, 73, 066301 

Bonanno A. & Urpin V. 2007, ApJ, 376, 214 

Campbell C.G. 1987, MNRAS, 229, 405 

Campbell C.G. 1992, GApFD, 63, 179 

Chandrasekhar S. 1960, Proc. Natl. Acad. Sci., 46, 253 

Davidson K. & Ostriker J. P 1973, ApJ, 179, 585 

Elstner D. & Riidiger G. 2000, A&A, 358, 612 

Ghosh P, Lamb F.K. & Pethick C.J. 1977, ApJ, 217, 578 

Ghosh P & Lamb F.K. 1979a, ApJ, 232, 259 

Ghosh P & Lamb FK. 1979b, ApJ, 234, 296 

Illarinov A. F. & Sunyaev R. A. 1975, Astr. Ap., 39, 185 

Kiziltan B. & Thorsett S.E. 2009, ApJ, 693, L109 

Kluzniak W. & Rappaport S. 2007, ApJ, 671, 1990 

Kulkarni A. K. & Romanova M. M. 2008, MNRAS, 386, 673 

Mestel L. 1961, MNRAS, 122, 473 

Miller K. A. & Stone J. M. 1997, 489, 890 

Romanova M. M., Ustyugova G. V., Koldoba A. V. & Lovelace R. V. E. 2002, 
ApJ, 578, 420 

Riidiger G., Hollerbach R., Schultz M. & Elstner D. 2007, MNRAS, 377, 1481 

Riidiger G. & Shalibkov D. 2002, A&A, 393, L81 

Shakura N.l. 1975, Soviet Astr. Letters, 1, 223 

Shakura N.l. & Sunyaev R.A. 1973, A&A, 24, 337 

Shalibkov D. & Riidiger G. 2000, MNRAS, 315, 762 

Tayler R.J. 1973, MNRAS, 161, 365 

Turner N.J. Sano T., 2008 ApJ, 679, L131 

Velikov E.P 1959, Sov. Phys. JETR 9, 995 



Wang Y.-M. 1987, A&A, 183, 257 



Appendix A: The Code 

In this Appendix we describe the computer code that we have 
used to solve equation ( |25] | and discuss some of the tests that we 
have performed to verify it. 

A.1. Description of tlie code 

In order to solve the 2-D elliptic PDE ( IZST i we use the Gauss- 
Seidel relaxation method. If we call the elliptic operator X 
and the right hand side b, then the original equation becomes: 
-C[S] = b. We turn this elliptic equation into a hyperbolic one by 
adding a pseudo time derivative; we can then consider the itera- 
tive procedure as a time evolution and write: dtS - -C[S] - b. In 
our case £^8^+ 45? - + ^) dy - ^ d, and b ^ 0. 

We then approximate the operators by discretizing the func- 
tions over a grid. The scheme that we use for discretizing the 
derivatives is as follows: 



deS\ij = 



drS\ij = 



2Aj 



and djSlij = 



2Ai 



and dJ.S\ij = 



■ 2Sij + 5, 



-(A.l) 



i-hj 



A/2 



(A.2) 



- S' 
At 



(A.3) 



where the indices / and j refer to grid points along the x and y 
coordinate directions respectively, while f represents the pseudo- 
time or iteration step. We use expressions ( IA.ll i-( IA.3l ) to dis- 
cretize equation dZSb and then, isolating the term >S'y , we get 
the iterative algorithm that we use in our code: 



1 S'.j^,-2S'^, + S'^} 



A/2 

'/+! 

I. J ■ 



A/ 



XiTI 



' cot vgrg\ S'l j^ i - S'. j_ j 



2Aj 



or C/+ 1 

v,-'-g'^,+i,;-'^,-i,; 



2Ai 



(A.4) 



We solve this proceeding outwards from / - I, j - 1; on the 
right hand side the terms that have already been calculated (i.e. 
the terms at positions / = / - 1 and j = y - 1) are taken at 
the cuiTent iteration t + 1, as usual in the Gauss-Seidel method. 
We provide an initial estimate for S at iteration zero and the 
code then modifies this by relaxing the solution using the chosen 
profiles of y,-, ve and rj, which are all functions of r and 6. 

The magnitude of the central dipole field and the accretion 
rate do not enter this equation directly, but they are used to cal- 
culate the inner edge of the disc ri„. We recall that we have a 
radial numerical domain going from to rout, and a physical 
domain of interest which begins at the same radius but stops ear- 
lier, at the light cylinder radius r\c . Along the direction the two 
domains coincide and go from 0top, which is the top surface of 
the corona, to 0eq, which is the equatorial plane (because of the 
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symmetry of the system with respect to the equatorial plane, the 
solution below the equator will be the same as above it). 

For the mass and radius of the neutron star, we use the canon- 
ical values, IAMq and 10 km respectively. We fix the accretion 
rate as m = 0.03 (in units of MEdd), giving a magnetospheric 
radius of about lOrg when Bq ~ 3 10^ G, as typical for a mil- 
lisecond pulsar. 

A.2. Testing of the code 

In order to check the code for stability and convergence, to esti- 
mate errors and to optimize the iteration procedure by choosing 
an appropriate iteration step, we performed a number of tests, 
some of which are now described. 

During this test phase we used the following values for the 
parameters: 

- magnetic field at the stellar surface: Bq = 3 10*^ G; 

- size of the domain: rin = 10 rg, rout = 750 rg, 0top - 80°, 
0eq = 90° and 0^ = 82°; 

- radial velocity at inner edge: Vr('"in) - vo = 10^ cm s"^ 
which is the value obtained from equation ( fTTI ) when 
(a, m) = (0.15, 0.03) or (0.1, 0.07); 

- diffusivity: r,,in = nn, r;,out = '"out, m = 10'" cm^ s"' and 
771; - 10'' cm^ s"'; 

- initial estimate for the magnetic stream function: 

>S = ro ■ sin 0/r°-5 (for a dipolar field >S'*'p = r,, sin^ 0/r); 

- iteration time step: Af = 3 10"^ ■ Ax ■ Ay. 

The tests can be divided into two main groups: with and 
without a known analytic solution. For the latter, we can estimate 
errors by calculating the residuals and comparing the solutions 
obtained with different grid resolutions, while for the former we 
also have the difference between the computed values and the 
exact results. 

A.2.1 . Test with an analytic solution 

There are two cases for which we can obtain an analytic solu- 
tion. The first has dipolar boundary conditions and no poloidal 
motion (v^ - vq - 0); in this case the poloidal component of 
the field must be dipolar everywhere (we refer to this test as D, 
for dipolar). In Section [3] we showed that, in order to have a 
dipolar field as a stationary solution of the induction equation, 
the poloidal velocity has to fulfill the relation given by equation 
(I2TI 1 and using v, - vq -Ois, consistent with that condition. Our 
second analytic test case has the boundary conditions for S set to 
zero. In this case, regardless of the profile used for the velocity, 
S(r, 0) = is a solution in all of the domain (we call this test 
Z, for zero). Even if at first glance a test with an identically zero 
solution may seem to be of little importance, we think that it is 
useful, because in this way we can test the code by including all 
of the terms that will be present when solving for the cases of 
interest (i.e. including v^, vq and rj). 

In both cases, we test two different configurations by chang- 
ing the velocity profile. In test D we consider a case with zero 
velocity (test Dl) and another one where vg is given by equation 
(l2n i and with given by equation ( fTTT l in the corona and being 
zero in the disc (test D2). In test Z we consider the same velocity 
profile as the one that we will use for our cases of interest, given 
by equations (fTTT l and ( |2TI ) (test Zl), and a velocity profile which 
is the same as that used in test D2 (test Z2). 

In all of these tests we follow the same procedure, we ver- 
ify the stability of the code, we estimate the error and see if it 
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Fig. A.l. Relative error at each grid point of the region of inter- 
est for a representative case: test D2 with a 400 x 20 grid. For 
the other grids (and tests) the results have exactly the same qual- 
itative behaviour, but the precise value of the error is different 
and increases if the number of grid-points is decreased. For the 
100 X 20 grid, the maximum relative error is ~ 15 per cent, while 
for the 200 x 20 grid it is about 0.8 per cent and finally for the 
400 X 20 grid (shown in the figure) it is 0.15 per cent. 



scales coiTectly, checking the convergence of the solution. We 
do this by studying how the numerical solution changes when 
varying the number of grid points (A^, and Nj) and the number of 
iterations. 

We use five grids in total. When testing the dependence on 
A^^ we use: 200 x 20, 200 x 40 and 200 x 80; while when testing 
the dependence on A^, we use: 100 x 20, 200 x 20 and 400 x 20. 
For each of these grids we calculate: (i) the absolute difference 
and (ii) the relative difference, between the numerical solution 
and the analytic one at each point of the grid; and (iii) the root 
mean square (rms) of the numerical solution S at each iteration 
step. The results obtained are very similar for all of the five grids 
and for each of the four tests and can be summarized with the 
following four statements: (1) both the absolute eiTor and the 
relative error have a maximum near to the inner edge and 
then decrease quite rapidly. With the 400 x 20 grid, the maxi- 
mum relative error is about 0. 1 per cent (see fig. lA. Il l; (2) chang- 
ing Nj does not produce any visible effect: while increasing Nj 
by a factor of 4 (from 20 to 80) decreases the maximum rela- 
tive error only slightly (~ 7 per cent), changing A^, from 100 
to 400 has a much greater effect, giving a decrease in the error 
of two orders of magnitude; (3) the reduction in the rms and in 
the maximum error becomes progressively smaller with increas- 
ing A^,, thus showing that we have convergence of the numerical 
solution; (4) using a sparser grid gives smaller errors at the be- 
ginning and during the relaxation process, however if one keeps 
iterating until the saturation level is reached, then the error with 
sparser grids is larger than with denser grids (suggesting that this 
problem could be suited for a multigrid approach). Regarding 
statements (2), (3) and (4), see fig. lA.2l 

A.2.2. Test with an unknown solution 

We use now a configuration with dipolar boundary conditions 
and a velocity field given by equations (fTTI) and (|2T]) . This is 
very similar to test D2, but in this configuration Vr is not taken 
to be zero in the disc. Even if we do not know the solution for 
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Fig.A.2. Maximum analytic error for the test D2 with three 
grids, which differ only in A^, . Increasing Nj from 20 to 80 pro- 
duces a very small decrease in the eiTor, only ~ 7%. 

this setup, we know from equation ( fTOl i that it has to approach 
a dipolar field when the coefficients — and — both go to zero. 
In order to test this we have considered five configurations, each 
with a different value of 770 ranging between 10" and 10'^ cm^ 
s"' (?7c is always taken to be two orders of magnitude larger than 
770). Figure lA3] shows clearly that for increasing 77, the rms of the 
numerical solution is approaching that for a dipole calculated on 
the same grid. 

For these five configurations we also performed the tests pre- 
viously described, i.e. the ones regarding changing the grid and 
comparing the errors and the rms. The results are again similar 
and confirm the four statements made earlier 

A.2.3. Other tests 

We next used the configuration of test D2 for checking three 
more aspects: determining the importance of the initial estimate 
for S and investigating which values to choose for the location 
of the outer radial boundary and for the iteration step. 

The kind of algorithm which we are using to solve equation 
([Tol l needs an initial estimate for the solution. According to how 
good or bad this estimate is with respect to the correct solution, 
one needs a smaller or larger number of iterations for completing 
the process. In order to show this and also to demonstrate that 
the final solution does not depend on the initial profile, we used 
four initial estimates for the magnetic stream function S: (1) a 
constant value, (2) a Gaussian profile (centred on r = 100 rg 
and with a width of 20 rg), (3) a profile increasing with r' (this 
gives B,- and Be increasing linearly with r) and (4) the profile 
which gives a dipolar magnetic field (the analytic solution for 
this configuration is - ro-sinS/r"-^). In all of the cases the final 
solution is the same, even for configuration (3), but the number 
of iterations required to reach saturation changes and goes from 
for case (4) to 2^" for cases (1) and (2) and to 2^^ for case (3). 

As mentioned previously (in Sections l2land lA. Il l, our region 
of physical interest goes from the inner edge of the disc rin out 
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Fig. A.3. The rms of the numerical solution for different values 
of 77 after 2^"^ iterations. When 77 ^ 00, the rms should reach the 
value for a dipole, which is plotted with a dotted line for both 
grids. 

to the light cylinder r\^. Since we do not want the solution in this 
region to be influenced by the outer boundary condition, we ran 
some tests using different values for the radius of the outer edge 
Tout and then compared the numerical solutions in the region of 
interest. We used the same setup as for the real problem that we 
were wanting to solve, i.e. with dipolar boundary conditions, the 
velocity field given by equations (fTTT l and (1211 1 and the diffusiv- 
ity given by equation ( fTSI l. We used six values of rout (150 r^, 
200 rg, 250 rg, 300 rg, 500 rg and 750 rg) and we found that the 
difference within the region of interest between the numerical 
solutions obtained using two subsequent values of rout became 
progressively smaller, until one could barely distinguish the dif- 
ferent solutions. We decided to put the outer boundary at 750 rg 
for the physical analysis; this gives results differing from those 
with rout — 500 rg by less than about 5 per cent. 

Finally we considered varying the iteration step size, i.e. the 
Af in equation (IA.4l i. that is written as c • Ar • A0. There is no 
simple argument of principle that can be used to determine the 
best value for c, therefore we determined it experimentally. We 
considered the same configuration as in test D2 and ran it several 
times varying only the value of c, going from 0.025 upwards. 
We found that the final asymptotic error was always the same, 
but that the number of iterations required to relax to the final 
solution was changing, decreasing as c increased. However there 
is an upper limit: when c > Cn,ax = 1 -25 the numerical solution 
diverges. Transferring this condition to Af, we obtain Af^^x = 
8.5 lO^^*. We can then change the way in which the iteration step 
size is calculated in the code and write: At - n ■ Af^ax, with n 
always smaller than 1. We find that using the value n - 0.95 is 
a good compromise in minimizing the number of iterations and 
preserving the code stability. 



